Note about prob, stat with R

stat_vs_prob

Chúng ta có thể sử dụng R base hoặc thư viện tidyverse để xử lý dữ liệu trong R.
Các ví dụ ở dưới sẽ cố gắng sử dụng cả 2 cách.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble  3.1.1     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ✓ purrr   0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
# Remove scientific notation
options(scipen=999)

Basic

cbind, rbind, apply

Giả sử ta có một bảng như thế này table

Có thể nhập dữ liệu theo dòng

stu1 <- c(73, 75, 74, 74)
stu2 <- c(95, 94, 12, 95)
stu3 <- c(66, 67, 63, 100)

# Sau khi nhập theo dòng, gộp các dòng này thành 1 dataframe
stu <- rbind(stu1, stu2, stu3) # rbind means combine rows
colnames(stu) <- c("Test1", "Test2", "Test3", "Test4") # rename column name
stu
##      Test1 Test2 Test3 Test4
## stu1    73    75    74    74
## stu2    95    94    12    95
## stu3    66    67    63   100

Có thể apply 1 function nào đó vào từng dòng/cột của dataframe

# Apply function mean, to each row of stu matrix.
means <- apply(stu, 1, mean) # 1 means row, 2 means column
medians <- apply(stu, 1, median)
stu <- cbind(stu, means, medians) # merge means and medians to original dataframe with cbind
stu
##      Test1 Test2 Test3 Test4 means medians
## stu1    73    75    74    74    74    74.0
## stu2    95    94    12    95    74    94.5
## stu3    66    67    63   100    74    66.5

quantile

x <- c(1, 4, 7, 9, 10, 14, 15, 16, 20, 21)
# Caluculate quantile 25%, 50% and 75% of vector
quantile(x, probs=c(0.25, 0.5, 0.75))
##   25%   50%   75% 
##  7.50 12.00 15.75

mean, variance

mean(x) # mean
## [1] 11.7
var(x) # variance (n-1)
## [1] 44.01111
sd(x) # standard deviation (n-1)
## [1] 6.634087
summary(x)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    7.50   12.00   11.70   15.75   21.00

correlation

Tìm correlation coefficient của 2 biến định lượng

# Using mtcars package to get data mtcars
head(mtcars, 5)
##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
# We check the relation between mpg and wt variables
mtcars %>%
  ggplot(aes(x = mpg, y = wt)) +
  geom_point()

Now, run cov to check if there is correlation between 2 variables

cor(mtcars$mpg, mtcars$wt)
## [1] -0.8676594

This indicates a negative linear relationship between the two variables.

Counting

Tổ hợp chập 3 của 5: chọn 3 phần tử trong 5 phần tử. Không quan tâm thứ tự

choose(5, 3)
## [1] 10

Sequence

Tạo ra 1 chuỗi từ 1 đến 100, với khoảng cách 4

seq(1, 100, 4)
##  [1]  1  5  9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93 97
c(1:100) # khoảng cách 2 phần tử là 1
##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##  [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
##  [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
##  [55]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
##  [73]  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
##  [91]  91  92  93  94  95  96  97  98  99 100

Set seed

# truyền vào 1 giá trị số bất kỳ, để xác định seed mình muốn tạo. Ví dụ ở đây là seed 1
set.seed(1)
# Ngay sau khi gọi set.seed thì gọi sample để tạo random sample cho seed 1 này
sample(LETTERS, 5)
## [1] "Y" "D" "G" "A" "B"

Bây giờ ta có thể mở rộng seed 1, bằng việc tạo sample với nhiều giá trị hơn

set.seed(1)
sample(LETTERS, 7)
## [1] "Y" "D" "G" "A" "B" "K" "N"

Dễ thấy, 5 giá trị đầu tiên của sample vẫn không thay đổi, nó chỉ bổ sung thêm 2 giá trị ngẫu nhiên vào seed này. Đây chính là ý nghĩa của việc dùng seed - đảm bảo khả năng reproduced của sample.

Bây giờ ta thử tạo 1 cái seed 2 với tập mẫu khác

set.seed(2)
sample(LETTERS, 7)
## [1] "U" "O" "F" "X" "H" "Q" "Z"

Sampling

Simple Random sampling

Lấy ngẫu nhiên n phần tử trong tập N. Xác suất để lấy mỗi phần tử là như nhau. Việc lấy mẫu ngẫu nhiên như này chỉ phù hợp nếu tổng thể có tính đồng nhất (homogeneous) vì khi đó 1 mẫu ngẫu nhiên mới có thể đại diện được cho tổng thể (phản ánh được 1 số thuộc tính của tổng thể).

Stratified Sampling

Chọn mẫu phân tổ là phương pháp lấy mẫu trong trường hợp tổng thể không đồng nhất. Cơ bản thì nếu tổng thể có thể chia nhỏ thành những tổng thể con, thì ta có thể sử dụng cách lấy mẫu này.
Ví dụ một số nhân tố có thể giúp chia nhóm mẫu là giới tính, khu vực địa lý, …

Systematic Sampling

Trường hợp ta biết tổng thể, và muốn tiết kiệm chi phí (tài nguyên) ta cần chọn mẫu từ tổng thể đó 1 cách có hệ thống (có quy luật).

Ví dụ tổng thể 100 phần tử, lấy mẫu 10 phần tử. Thì quy luật là trong 10 phần tử đầu tiên, random 1 cái. Sau đó cứ cộng 10. Ví dụ random 7, thì sẽ lấy mẫu 7, 17, 27, …, 97

Cluster Sampling

Toàn bộ tập quần thể sẽ được chia thành từ cụm hoặc thành từng phần. Sau đó chúng ta sẽ chọn ngẫu nhiên từng cụm. Tất cả các cá thể trong cụm đó sẽ được sử dụng làm tập mẫu.

The main difference between cluster sampling and stratified sampling is that in cluster sampling, the cluster is treated as the sampling unit and analysis is done on a population of clusters

Exploring Data

Frequency table

Dữ liệu ví dụ: data về vé phạt tại Philadelphia.

tickets <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-12-03/tickets.csv")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   violation_desc = col_character(),
##   issue_datetime = col_datetime(format = ""),
##   fine = col_double(),
##   issuing_agency = col_character(),
##   lat = col_double(),
##   lon = col_double(),
##   zip_code = col_double()
## )
summary(tickets)
##  violation_desc     issue_datetime                     fine        
##  Length:1260891     Min.   :2017-01-01 00:00:00   Min.   :  15.00  
##  Class :character   1st Qu.:2017-04-05 00:37:30   1st Qu.:  26.00  
##  Mode  :character   Median :2017-07-05 09:11:00   Median :  36.00  
##                     Mean   :2017-07-03 06:07:59   Mean   :  45.41  
##                     3rd Qu.:2017-10-01 13:20:30   3rd Qu.:  51.00  
##                     Max.   :2017-12-31 15:42:00   Max.   :1001.00  
##                                                                    
##  issuing_agency          lat             lon            zip_code     
##  Length:1260891     Min.   :39.57   Min.   :-75.99   Min.   :19102   
##  Class :character   1st Qu.:39.95   1st Qu.:-75.18   1st Qu.:19106   
##  Mode  :character   Median :39.95   Median :-75.16   Median :19123   
##                     Mean   :39.97   Mean   :-75.16   Mean   :19124   
##                     3rd Qu.:39.97   3rd Qu.:-75.15   3rd Qu.:19143   
##                     Max.   :40.37   Max.   :-74.96   Max.   :19154   
##                                                      NA's   :173588

Ta sẽ phân tích một số cột nhất đinh, là violation_desc, issue_datetime, fine và issuing_agency

tickets <- tickets %>%
  select(violation_desc, issue_datetime, fine, issuing_agency)
tickets %>%
  head(5)
## # A tibble: 5 x 4
##   violation_desc      issue_datetime       fine issuing_agency
##   <chr>               <dttm>              <dbl> <chr>         
## 1 BUS ONLY ZONE       2017-12-06 12:29:00    51 PPA           
## 2 STOPPING PROHIBITED 2017-10-16 18:03:00    51 PPA           
## 3 OVER TIME LIMIT     2017-11-02 22:09:00    26 PPA           
## 4 OVER TIME LIMIT     2017-11-05 20:19:00    26 PPA           
## 5 STOP PROHIBITED CC  2017-10-17 06:58:00    76 PPA

Đầu tiên hãy bắt đầu với câu hỏi, những lý do vi phạm violation_desc là những lý do gì và tần suất các lý do. Cái nào phổ biến nhất, ít phổ biến nhất.

# R base | %>% is from tidyverse, we use it here just to display the result
xtabs(~tickets$violation_desc) %>% 
  head(5)
## tickets$violation_desc
## +4HR IN LOADING ZONE    BLOCKING DRIVEWAY BLOCKING DRIVEWAY CC 
##                    1                 7546                  958 
## BLOCKNG MASS TRANSIT BUS NOT IN BUS STAND 
##                  182                   66
# tidyverse
tickets %>%
  group_by(violation_desc) %>%
  count() %>%
  head(5)
## # A tibble: 5 x 2
## # Groups:   violation_desc [5]
##   violation_desc           n
##   <chr>                <int>
## 1 +4HR IN LOADING ZONE     1
## 2 BLOCKING DRIVEWAY     7546
## 3 BLOCKING DRIVEWAY CC   958
## 4 BLOCKNG MASS TRANSIT   182
## 5 BUS NOT IN BUS STAND    66

Để biết tỷ trọng của từng loại vi phạm (chiếm bao nhiêu % trong tổng số vi phạm) thì sử dụng prop.table

# R base
prop.table(
  xtabs(~tickets$violation_desc)
) %>%
  `*`(100) %>% 
  round(2) %>%
  head(5)
## tickets$violation_desc
## +4HR IN LOADING ZONE    BLOCKING DRIVEWAY BLOCKING DRIVEWAY CC 
##                 0.00                 0.60                 0.08 
## BLOCKNG MASS TRANSIT BUS NOT IN BUS STAND 
##                 0.01                 0.01
# tidyverse
tickets %>%
  group_by(violation_desc) %>%
  summarise(freq = n()) %>%
  mutate(prop = round(freq / sum(freq) * 100, 4)) %>%
  # chỉ lấy tỷ lệ lớn hơn 1%
  filter(prop > 1)
## # A tibble: 15 x 3
##    violation_desc         freq  prop
##    <chr>                 <int> <dbl>
##  1 BUS ONLY ZONE         19867  1.58
##  2 EXPIRED INSPECTION   138575 11.0 
##  3 FIRE HYDRANT          20945  1.66
##  4 HP RESERVED SPACE     14554  1.15
##  5 LOADING ZONE   CC     12826  1.02
##  6 METER EXPIRED        181329 14.4 
##  7 METER EXPIRED CC     281060 22.3 
##  8 OVER TIME LIMIT      156859 12.4 
##  9 OVER TIME LIMIT CC    24585  1.95
## 10 PARKING PROHBITED     47232  3.75
## 11 PARKING PROHBITED CC  45082  3.58
## 12 PASSENGR LOADNG ZONE  24359  1.93
## 13 SIDEWALK              18716  1.48
## 14 STOP PROHIBITED CC   115898  9.19
## 15 STOPPING PROHIBITED   47395  3.76

Vẽ đồ thị - Bar plot

# R base | not recommend for plotting
prop.table(
  xtabs(~tickets$violation_desc)
) %>%
  `*`(100) %>% 
  round(2) %>%
  barplot()

# tidyverse
tickets %>%
  group_by(violation_desc) %>%
  summarise(freq = n()) %>%
  mutate(prop = round(freq / sum(freq) * 100, 4)) %>%
  # chỉ lấy tỷ lệ lớn hơn 1%
  filter(prop > 1) %>%
  ggplot(aes(x = reorder(violation_desc, freq), y = freq)) +
    geom_bar(stat="identity") +
    coord_flip()

Frequency table with 2 category variables

Vẫn sử dụng data vilation ở trên

xtabs(~tickets$violation_desc + tickets$issuing_agency) %>% 
  head(5)
##                       tickets$issuing_agency
## tickets$violation_desc CENTER C FAIRMNT HOUSING PENN POLICE POST OFF  PPA
##   +4HR IN LOADING ZONE        0       0       0    0      1        0    0
##   BLOCKING DRIVEWAY           0       3       4   16   6066        5 1369
##   BLOCKING DRIVEWAY CC        1       0       0   22    215        0  720
##   BLOCKNG MASS TRANSIT       30       0       0    3     91        0   20
##   BUS NOT IN BUS STAND        0       0       0    0      4        0   62
##                       tickets$issuing_agency
## tickets$violation_desc PRISONS SEPTA TEMPLE
##   +4HR IN LOADING ZONE       0     0      0
##   BLOCKING DRIVEWAY         47     1     35
##   BLOCKING DRIVEWAY CC       0     0      0
##   BLOCKNG MASS TRANSIT       0    37      1
##   BUS NOT IN BUS STAND       0     0      0

Distribution

Binom Distribution

Phân bố binomial thường được dùng để mô tả những biến ngẫu nhiên nhị phân (xảy ra vs không xảy ra)

Nếu ta thực hiện n phép thử nhị phân, với xác suất xảy ra True là p. Thì số lần xảy ra True kỳ vọng là $np

# Ví dụ sample có n = 20 phép thử, pi là 0.5 là xác suất để sự kiện xảy ra.
# Thực hiện 10000 lần cái sample n kia và tính xác suất
bino.gen(samples = 10000, n = 20, pi = 0.5)

## $simulated.distribution
## Successes
##     2     3     4     5     6     7     8     9    10    11    12    13    14 
## 0.000 0.001 0.005 0.014 0.038 0.074 0.123 0.161 0.179 0.162 0.112 0.078 0.032 
##    15    16    17    18    19 
## 0.016 0.004 0.001 0.000 0.000 
## 
## $theoretical.distribution
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 0.000 0.000 0.000 0.001 0.005 0.015 0.037 0.074 0.120 0.160 0.176 0.160 0.120 
##    13    14    15    16    17    18    19    20 
## 0.074 0.037 0.015 0.005 0.001 0.000 0.000 0.000

Có thể thấy, nếu thực hiện rất nhiều sample, thì phân bố binomial sẽ tiệm cận normal với giá trị kỳ vọng là \(np\)

Poisson Distribution

Phân phối Poisson thường được dùng để mô tả những biến ngẫu nhiên rời rạc, dạng biến đếm.

Ví dụ tần suất 1 sự kiện xảy ra trong 1 khoảng thời gian.

  • Số lượng xe ô tô đi qua 1 điểm trong khoảng thời gian xác định (ví dụ 30s)

  • Số lượng khách hàng vào cửa hàng trong 1 khoảng thời gian xác định (ví dụ từ 9h đến 10h sáng)

Như vậy, biến ngẫu nhiên X phân phối Poisson với tham số \(\lambda\), ở đây \(\lambda\) là giá trị trung bình.

\[ X \sim {\sf Poisson}(\lambda) \]

Ví dụ, biến ngẫu nhiên X, có phân phối Poisson và \(\lambda\) là 3, tức là giá trị trung bình của nó là 3.

Vẽ biểu đồ ra, thì điểm cao nhất của phân bố sẽ xoay quanh 3

hist(rpois(100000, 5), breaks = 15, col = 'orange2')

Có thể tính CDF, \(P(X < \alpha)\) bằng cách

# Với lambda là 5, muốn tính xác suất X nhỏ hơn 8 ta dùng
ppois(8, 5)
## [1] 0.9319064

Normal distribution

normal

Chi-Square Distributions

Khi ta có một biến ngẫu nhiên có phân bố chuẩn tắc là Z với \(\mu = 0\) và \(\sigma = 1\), khi đó bình phương của Z sẽ có phân bố chi-square

chi-square Nếu chỉ có 1 biến là \(Z^2\) thì phân bố Chi-Squared sẽ có 1 bậc tự do.

chi-square-df

Nếu có tổng nhiều \(Z^2\) thì số bậc tự do cũng tăng lên (bằng số biến \(Z^2\))

chi-square-more ## Student’s t Distributions

Kết hợp 1 biến phân bố chuẩn tắc, và 1 biến phân bố Chi-squared ta được 1 biến phân bố student như sau

student

F Distributions

2 biến phân bố Chi-squared sẽ có tỷ lệ là 1 phân bố F.

f

Statistics from sample - Đặc trưng / tham số của mẫu

Đặc trưng từ mẫu phân phối chuẩn

Lấy mẫu ngẫu nhiên \(n\) phần tử trong tổng thể có phân phối Chuẩn.

Mỗi quan sát trong mẫu này là 1 biến ngẫu nhiên \(X_i\) có phân phối Chuẩn với trung bình là \(\mu\) và phương sai \(\sigma^2\)

stat_normal

Đặc trưng từ mẫu phân phối Poisson

poisson

Point estimation

Ước lượng điểm. Ví dụ với 1 biến ngẫu nhiên phân phối chuẩn với tham số \(\mu\) và \(\sigma\), ta lấy mẫu ngẫu nhiên, và từ mẫu ngẫu nhiên N phần tử này ta cố gắng ước lượng 2 tham số ở trên.

Có thể đơn giản là lấy kỳ vọng mẫu làm ước lượng cho \(\mu\) tổng thể và phương sai mẫu làm ước lượng cho \(\sigma\) tổng thể.

Giả định, ta có 1 ước lượng điểm \(T\) của 1 parameter \(\theta\)

Các tính chất của Point estimation

Mean Square Error

Với 1 phương pháp lấy ước lượng điểm, thực hiện trên nhiều sample thì ta có các ước lượng \(T_1\), \(T_2\), \(T_3\), …

Đây là 1 chỉ số để đánh giá độ chuẩn xác của 1 ước lượng điểm. Nó là bình phương của sai số giữa \(T\) và \(\theta\) và sau đó thì tính mean.

\[ \begin{aligned} MSE[T] &= E[(T - \theta)^2] \\ &= Var[T] + (E[T] - \theta)^2 \\ &= Var[T] + (Bias(T))^2 \end{aligned} \] MSE càng nhỏ thì ước lượng càng chính xác.

Unbiased Estimators

Có thể thấy 1 trong 2 thành phần của MSE là Bias của T chính là hiệu của \(E[T]\) và \(\theta\)

Unbiased Estimators là estimator T mà \(E[T] = \theta\)

Sample mean \(\overline{X}\) và sample variance \(S^2\) là unbiased estimators của population mean và population variance

Efficiency

Một estimator tốt là 1 estimator có MSE nhỏ, và Efficiency cũng là 1 cách để nói về MSE nhỏ

\[ eff[T] = \frac{1}{MSE[T]} \]

Consistency

Với các sample size khác nhau, ta có thể tạo ra các estimator khác nhau. Ví dụ, biến ngẫu nhiên \(N(\mu, \sigma)\) lấy sample N = 1.

Khi đó, ước lượng \(T_1\) được tính từ sample 1 phần tử trên.

Thêm 1 quan sát vào sample trên, ta có sample 2 phần tử, từ đó tính được \(T_2\), tương tự thế, tăng số lượng quan sát, ta có các \(T_n\) ước lượng, đối với sample có n phần tử/quan sát.

Ước lượng là consistent nếu n càng tăng, thì ước lượng \(T\) càng tiến đến/hội tụ về giá trị thực sự \(\theta\)

Robustness

Một chỉ số phản ánh mức độ phụ thuộc của estimator vào phân phối của tổng thể.

Nói cách khác, trường hợp mà ta không biết tổng thể có phân phối gì, thì việc estimate sẽ khó và thiếu chính xác hơn. Nên là 1 estimate mà càng ít phụ thuộc vào loại phân phối của tổng thể thì càng tốt, càng Robust.

Mở rộng, thì ngoài phân phối, các giả thiết khác như tính độc lập, … cũng được xét đến trong Robustness

Cách xác định Point estimator

Có 2 phương pháp chính dùng để tính ước lượng điểm.

  • Method of moments: Dựa vào hàm tạo moment (Moment generate function) của phân phối.
  • MLE: Maximum likelyhood estimate

Dựa vào Method of moments mà ta xác định được 1 estimator của population mean là sample mean. (vì sample mean là 1st moment)

Ví dụ về MLE:

MLE cho Normal distribution MLE cho Poisson distribution

Fisher Information

Thông tin Fisher - sẽ tìm hiểu thêm sau khi cần thiết.

Confidence Interval

Khoảng tin cậy (CI) quan trọng ở cách tính.

Với 1 cách tính khoảng tin cậy 95% xác định, với 1 sample X xác định thì ta tính được khoảng tin cậy \(L_x\) đến \(U_x\); có nghĩa là nếu ta lấy thêm 100 samples nữa, Y, Z, K, … thì 95 samples trong số đó, khi áp dụng cách tính khoảng tin cậy trên, ra được các khoảng \(L_y, U_y\); \(L_k, U_k\) chứa true parameter.

CI for Population Means

Sampling từ phân phối chuẩn và biết \(\sigma^2\) của tổng thể

\[ \left[ \overline{x} - z_{1 - \alpha/2}\frac{\sigma}{\sqrt{n}}, \; \overline{x} + z_{1 - \alpha/2}\frac{\sigma}{\sqrt{n}} \right] \] Có thể vẽ CI đối với trường hợp trên như sau,

set.seed(10)
cisim(samples = 100, n = 100, parameter = 0.5, type = "Mean", sigma = 1, conf.level = 0.95)

## 5 % of the random confidence intervals do not contain Mu = 0.5 .

Trong đó, \(z_1-\alpha/2\) có thể được tính bằng

# với alpha là 5% thì giá trị z là
qnorm(1 - 0.05/2)
## [1] 1.959964

Để tính CI, có thể áp dụng công thức trên, hoặc dùng z.test

z.test(
  x = c(100, 120, 130, 140, 150, 145, 134, 123, 156, 172), # Data từ sample
  sigma.x = 3, # sigma của tổng thể
  n.x = 10, # size của sample
  conf.level = .95 # độ tin cậy 95%
)
## 
##  One Sample z-test
## 
## data:  c(100, 120, 130, 140, 150, 145, 134, 123, 156, 172)
## z = 144.41, p-value < 0.00000000000000022
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  135.1406 138.8594
## sample estimates:
## mean of x 
##       137

Sampling từ phân phối chuẩn và không biết \(\sigma^2\) của tổng thể

Thay vì sử dụng \(\sigma\) tổng thể, ta dùng \(s\) là sample standard deviation

\[ \left[ \overline{x} - t_{1 - \alpha/2;n-1}\frac{s}{\sqrt{n}}, \; \overline{x} + t_{1 - \alpha/2;n-1}\frac{s}{\sqrt{n}} \right] \]